#include "Calibration.hh"
#include "TMath.h"

using namespace std;

Calibration::Calibration(){
};
Calibration::Calibration(const char* filename, Settings* setting){
  SetFile( filename );
  fSett = setting;
  fVerbose = fSett->VLevel();
  fNofDgfs = fSett->NofDgfs();
  fNofDgfChans = fSett->NofDgfChans();
  fNofAdcs = fSett->NofAdcs();
  fNofAdcChans = fSett->NofAdcChans();
  fBeamdumpDgf = fSett->BeamdumpDgf();
  ReadCalibration();
  fRand = new TRandom();
  if(fVerbose>0)
    PrintCalibration();
};
Calibration::~Calibration(){
  cout << "destructor" << endl;
};
void Calibration::ReadCalibration(){

  TEnv *config = new TEnv(fInputFile.data());
  if(fVerbose>1)
    cout << "reading dgfs" << endl;
  fDgfOffset.resize(fNofDgfs);
  fDgfGain.resize(fNofDgfs);
  for(int dgf=0; dgf<fNofDgfs; dgf++){
    fDgfOffset[dgf].resize(fNofDgfChans);
    fDgfGain[dgf].resize(fNofDgfChans);
    for(int chan=0; chan <fNofDgfChans; chan++){
      fDgfOffset[dgf][chan] = config->GetValue(Form("dgf_%d_%d.Offset", dgf, chan),0.);
      fDgfGain[dgf][chan] = config->GetValue(Form("dgf_%d_%d.Gain", dgf, chan),1.);
    }
  }
  if(fVerbose>1)
    cout << "reading beamdump" << endl;
  fBeamdumpOffset = config->GetValue(Form("dgf_%d_%d.Offset", fBeamdumpDgf, 0),0.);
  fBeamdumpGain = config->GetValue(Form("dgf_%d_%d.Gain", fBeamdumpDgf, 0),1.);

  if(fVerbose>1)
    cout << "reading adcs" << endl;
  fAdcOffset.resize(fNofAdcs);
  fAdcGain.resize(fNofAdcs);
  for(int adc=0; adc<fNofAdcs; adc++){
    fAdcOffset[adc].resize(fNofAdcChans);
    fAdcGain[adc].resize(fNofAdcChans);
    for(int chan=0; chan<fNofAdcChans; chan++){
		fAdcOffset[adc][chan] = config->GetValue(Form("adc_%d_%d.Offset", adc, chan),0.);
      fAdcGain[adc][chan] = config->GetValue(Form("adc_%d_%d.Gain", adc, chan),1.);      
    }
  }
  if(fVerbose>1)
    cout << "reading cd" << endl;

  for(int quadr=0; quadr<4; quadr++){
    for(int strip=0; strip<16; strip++){
      fCDRingGain[quadr][strip] = config->GetValue(Form("CD.Quadrant%d.ring.%d.gain",quadr,strip),1.);
      fCDRingOffset[quadr][strip] = config->GetValue(Form("CD.Quadrant%d.ring.%d.offset",quadr,strip),0.);
      fCDStripGain[quadr][strip] = config->GetValue(Form("CD.Quadrant%d.strip.%d.gain",quadr,strip),1.);
      fCDStripOffset[quadr][strip] = config->GetValue(Form("CD.Quadrant%d.strip.%d.offset",quadr,strip),0.);
    }
  }
  if(fVerbose>1)
    cout << "reading tilt correction" << endl;
  for(int quadr=0; quadr<4; quadr++){
    fFTilt[quadr] = config->GetValue(Form("Forward.TiltCorr.%d",quadr),0.);
    fBTilt[quadr] = config->GetValue(Form("Backward.TiltCorr.%d",quadr),0.);
    fFTiltOffset[quadr] = config->GetValue(Form("Forward.TiltOffset.%d",quadr),0.);
    fBTiltOffset[quadr] = config->GetValue(Form("Backward.TiltOffset.%d",quadr),0.);
  }  
delete config;
};
void Calibration::PrintCalibration(){
  cout << "DGFs" << endl;
  for(int dgf=0; dgf<fNofDgfs; dgf++){
    for(int chan=0; chan<fNofDgfChans; chan++){
      cout << Form("dgf_%d_%d.Offset\t", dgf, chan) << fDgfOffset[dgf][chan] << endl;
      cout << Form("dgf_%d_%d.Gain\t", dgf, chan) << fDgfGain[dgf][chan] << endl;
    }
  }
  cout << "Beamdump" << endl;
  cout << Form("dgf_%d_%d.Offset\t", fBeamdumpDgf, 0) << fBeamdumpOffset << endl;
  cout << Form("dgf_%d_%d.Gain\t", fBeamdumpDgf, 0) << fBeamdumpGain << endl;
  cout << "ADCs" << endl;
  for(int adc=0; adc<fNofAdcs; adc++){
    for(int chan=0; chan<fNofAdcChans; chan++){
      cout << Form("adc_%d_%d.Offset\t", adc, chan) << fAdcOffset[adc][chan] << endl;
      cout << Form("adc_%d_%d.Gain\t", adc, chan) << fAdcGain[adc][chan] << endl;
    }
  }
  cout << "CD" << endl;
  for(int quadr=0; quadr<4; quadr++){
    for(int strip=0; strip<16; strip++){
      cout << Form("CD.Quadrant%d.ring.%d.gain\t",quadr,strip) << fCDRingGain[quadr][strip] <<endl;
      cout << Form("CD.Quadrant%d.ring.%d.offset\t",quadr,strip) << fCDRingOffset[quadr][strip] <<endl;
      cout << Form("CD.Quadrant%d.strip.%d.gain\t",quadr,strip) << fCDStripGain[quadr][strip] <<endl;
      cout << Form("CD.Quadrant%d.strip.%d.offset\t",quadr,strip) << fCDStripOffset[quadr][strip] <<endl;
    }
  }
  cout << "Tilt corr" << endl;
  for(int quadr=0; quadr<4; quadr++){
    cout << Form("Forward.TiltCorr.%d\t",quadr) << fFTilt[quadr] << endl;
    cout << Form("Backward.TiltCorr.%d\t",quadr) << fBTilt[quadr] << endl;
    cout << Form("Forward.TiltOffset.%d\t",quadr) << fFTiltOffset[quadr] << endl;
    cout << Form("Backward.TiltOffset.%d\t",quadr) << fBTiltOffset[quadr] << endl;  }
};

double Calibration::DgfEnergy(int dgf, int chan, unsigned short raw){
  if( (dgf>-1) && (dgf<fNofDgfs) && (chan>-1) && (chan<fNofDgfChans) ){
    //cout << "dgf: "<<dgf<<" channel: "<<chan<<" gain: "<<setprecision(5)<<fDgfGain[dgf][chan]<<" offset: "<<fDgfOffset[dgf][chan]<<endl;
    return fDgfGain[dgf][chan]*( raw +0.5 - fRand->Uniform() - fDgfOffset[dgf][chan]);
  }
  else{
    cerr << "dgf "<<dgf<<" channel "<<chan<<" not found!" << endl;
    return -1;
  }
};

double Calibration::RemoveOffset(int forward, int quadr, unsigned short raw){
  if(forward==1){
    return  raw+0.5 - fRand->Uniform() - fFTiltOffset[quadr];
  }
  else if(forward==0){
    return  raw+0.5 - fRand->Uniform() - fBTiltOffset[quadr];
  }
  else{
    cerr << "forward (0,1) quadr (0,1,2,3) not correct: "<<forward<<"  "<<quadr<< endl;
    return -1; 
  }
};
double Calibration::DeltaEnergy(int adc, int chan, double rear){
  if( (adc>-1) && (adc<fNofAdcs) && (chan>-1) && (chan<fNofAdcChans) ){
    return fAdcGain[adc][chan]*(rear- fAdcOffset[adc][chan]);
  }
  else{
    cerr << "adc "<<adc<<" channel "<<chan<<" not found!" << endl;
    return -1; 
  }
};
double Calibration::AdcEnergy(int adc, int chan, unsigned short raw){
  if( (adc>-1) && (adc<fNofAdcs) && (chan>-1) && (chan<fNofAdcChans) ){
    return fAdcGain[adc][chan]*( raw +0.5 - fRand->Uniform() - fAdcOffset[adc][chan]);
  }
  else{
    cerr << "adc "<<adc<<" channel "<<chan<<" not found!" << endl;
    return -1; 
  }
};
/*
double Calibration::StripPos(int stradc, int strchan, unsigned short strraw, int rearadc, int rearchan, unsigned short rearraw){
  if( (stradc>-1) && (stradc<fNofAdcs) && (strchan>-1) && (strchan<fNofAdcChans) && (rearadc>-1) && (rearadc<fNofAdcs) && (rearchan>-1) && (rearchan<fNofAdcChans)){
    return fAdcGain[stradc][strchan]*( (double) strraw - fAdcOffset[stradc][strchan] )/( (double) rearraw - fAdcOffset[rearadc][rearchan] );
  }
  else{
    cerr << "strip adc "<<stradc<<" channel "<<strchan<<" not found!" << endl;
    cerr << "rear adc "<<rearadc<<" channel "<<rearchan<<" not found!" << endl;
    return -1; 
  }
};
*/
double Calibration::StripPos(int stradc, int strchan, unsigned short strraw, double rear){
  if( (stradc>-1) && (stradc<fNofAdcs) && (strchan>-1) && (strchan<fNofAdcChans) && (rear>0)){
	 return fAdcGain[stradc][strchan]*( (double) strraw - fAdcOffset[stradc][strchan] )/rear;
//	 return fAdcGain[stradc][strchan]*( (double) strraw - fAdcOffset[stradc][strchan] );
  }
  else{
//    cerr << "strip adc "<<stradc<<" channel "<<strchan<<" not found!" << " or rear < 0 " << rear << endl;
    return -1; 
  }
};
double Calibration::RearCorrection(double rear, double pos, int quadr, bool forward){
  double corrected;
  if(forward){
    corrected = rear/(fFTilt[quadr]*(1-pos)+1);
    if(fVerbose>0)
      if(rear>600 && rear<1200 )
	cout << "correcting " << rear << " -> " << corrected << " (pos) " << pos << " Tilt " << fFTilt[quadr] << " correction = " << corrected -rear << endl;
  }
  else{
    corrected = rear/(fBTilt[quadr]*(1-pos)+1);
	 if(fVerbose>0)
      if(rear>0 && rear<3300 )
	cout << "correcting " << rear << " -> " << corrected << " (pos) " << pos << " Tilt " << fBTilt[quadr] << " correction = " << corrected -rear << endl;
  }
  return corrected;
};

double Calibration::CDRingEnergy(int quadr, int ring, unsigned short raw){
  if( (quadr>-1) && (quadr<4) && (ring>-1) && (ring<16) ){
    if(quadr==0 && ring==5 && fVerbose>3)
      cout << raw << " -> " << fCDRingGain[quadr][ring]*( raw +0.5 - fRand->Uniform() - fCDRingOffset[quadr][ring]) << endl;
    return fCDRingGain[quadr][ring]*( raw +0.5 - fRand->Uniform() - fCDRingOffset[quadr][ring]);
  }
  else{

    cerr << "quadr "<<quadr<<" ring "<<ring<<" not found!" << endl;
    return -1; 
  }
};
double Calibration::CDStripEnergy(int quadr, int strip, unsigned short raw){
  if( (quadr>-1) && (quadr<4) && (strip>-1) && (strip<16) ){
    return fCDStripGain[quadr][strip]*( raw +0.5 - fRand->Uniform() - fCDStripOffset[quadr][strip]);
  }
  else{

    cerr << "quadr "<<quadr<<" strip "<<strip<<" not found!" << endl;
    return -1; 
  }
};
